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I.   INTRODUCTION 

In  October  of  1968,  Tropical  Storm  Kara  forced  R/V  PILLSBURY 
south  from  an  operating  area  in  the  vicinity  of  latitude  30°25'N, 
longitude  76°52'W.   While  in  the  stages  of  maturation,  the  tropical 
storm  displayed  the  curious  behavior  of  slowing  and  crossing  its  own 
path  several  times.   Kara  eventually  reached  hurricane  force,  increased 
speed  and  moved  on  a  generally  northeast  course.   In  transit  to  the 
operating  area  after  passage  of  the  storm,  the  research  vessel 
experienced  set  to  the  right  and  then  left,  indicating  that  a  hurricane- 
induced  eddy  had  been  crossed.   Bathythermograph  stations  were  taken 
in  order  to  locate  the  associated  dome-like  density  structure,  and  STD 
stations  across  the  presumed  dome  were  taken  after  it  had  been  located. 
The  track  of  Kara  and  the  track  of  R/V  PILLSBURY  are  shown  in  Figure  1. 
The  hurricane  track  was  obtained  from  the  advisories  issued  by  the 
National  Hurricane  Center. 

Kara  may  be  considered  to  have  been  quasi-stationary  during  the 
period  14-15  October  1968.   Therefore  it  affords  an  opportunity  to 
study  some  aspects  of  hurricane-induced  dynamics  through  the  upper 
part  of  the  water  column.   As  Stevenson  (1966)  points  out,  hurricanes 
are  unique  "laboratories"  for  the  study  of  air-sea  interaction  and 
ocean  dynamics.   Although  the  significant  atmospheric  effects  of  a 
hurricane  are  spread  over  horizontal  scales  on  the  order  of  several 
hundred  kilometers,  the  most  interesting  oceanic  effects  take  place 
within  an  area  which  can  be  crossed  in  one  day.   Unfortunately, 
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investigations  of  the  response  to  hurricanes  depend  by  necessity  on 
"cruises  of  opportunity."  Given  the  availability  of  research  vessels, 
studies  of  this  sort  are  further  complicated  by  the  unpredictability 
of  the  occurrence,  duration  and  track  of  a  storm. 

A  study  of  low  intensity  storms  such  as  Kara  is  particularly 
relevant  because  the  wind  velocities  are  not  too  different  from  those 
found  in  the  winter  storms  in  the  North  Atlantic  Ocean.   In  particular, 
the  data  provide  a  means  of  testing  assumptions  of  conservation  of 
potential  vorticity  in  the  ocean.  An  analysis  of  those  areas  where 
the  conservation  equations  are  not  applicable  is  important  because 
understanding  the  causes  of  the  differences  between  observed  and 
measured  response  contributes  to  an  understanding  of  the  response  of 
the  ocean  in  the  near  surface  layer.  A  third  reason  for  studying  such 
storms  is  that  the  deep  ocean  response  may  be  typical  of  numerous 
storms  of  near  hurricane  force.   If  so,  the  associated  changes  in  mass 
structure  are  significant  at  great  depths.   For  example,  the 
perturbation  of  the  isopycnals  will  have  drastic  effects  on  acoustic 
propagation  at  depths  at  least  as  great  as  the  SOFAR  channel.   The 
analysis  also  serves  as  an  addition  to  a  catalogue  of  techniques  useful 
in  analyzing  transient  oceanic  behavior. 

Circular  symmetry  is  assumed  in  the  model  used  in  this  paper. 
However,  it  can  be  seen  from  Figure  1  that  insufficient  stations  have 
been  taken  to  positively  identify  the  center  of  the  vortex.   Station  24 
was  made  after  the  vortex  had  been  crossed  once,  when  the  research 
vessel  returned  to  take  a  center  station.   It  was  felt  that  the  station 
was  taken  close  to  the  center  but  the  possibility  that  it  was  not 


cannot  be  excluded.   The  stationarity  of  Kara  increases  the  amplitude 
of  the  response.   Since  the  data  were  collected  several  days  after  the 
storm  passed, it  may  be  assumed  that  most  of  the  energy  is  found  in 
geostrophic  motion. 


II.   METHODOLOGY 

The  existence  of  a  cold  water  dome  after  the  passage  of  a  hurricane 
is  suggestive  of  upwelling.   This  was  concluded  by  Leipper  (1967)  and 
is  supported  by  data  collected  by  Landis  (1966).   Leipper's  data 
inspired  O'Brien  and  Reid  (1967)  to  pursue  their  study  of  upwelling  in 
a  two  layer  ocean  as  a  response  to  a  stationary  hurricane.  O'Brien 
and  Reid  established  through  a  transient,  non-linear  model  that  the 
wind  stress  associated  with  storms  of  hurricane  force  and  structure  is 
capable  of  producing  velocity  divergence  in  the  surface  layer  and  large 
scale  upwelling.   O'Brien  (1967)  added  turbulent  mixing  to  the  model 
and  produced  results  which  bear  resemblance  to  Leipper's  observations 
of  the  near  surface  layer. 

The  problem  can  be  approached  somewhat  differently.   It  is  well 
known  that  surface  divergence  and  associated  upwelling  will  stretch 
water  columns.  A  fluid  will  respond  by  spinning,  thereby  generating 
relative  vorticity  about  the  vertical  axis.   The  effect  extends  to 
the  bottom  of  the  fluid  and  is  transmitted  much  more  efficiently  than 
momentum  transfer  by  friction.   The  model  in  this  study  utilizes  the 
conservation  of  potential  vorticity  to  describe  the  interior  response. 
The  assumption  of  a  frictionless  interior  makes  the  continuously 
stratified  problem  tractable  by  allowing  use  of  the  conservation 
equation  in  the  interior.   This  model  and  O'Brien's  differ  because 
they  investigate  different  aspects  of  the  response  to  axisymmetric 
forcing.   O'Brien  solves  the  governing  equations  numerically  as  an 
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initial  value  problem  in  time  because  he  is  interested  in  the  evolution 
of  the  response.   By  its  very  nature,  his  two  layer  model  is  incapable 
of  obtaining  the  structure  of  the  interior.   However,  it  is  successful 
in  duplicating  many  aspects  of  the  time  dependent  response,   This 
study  seeks  to  model  the  steady-state  response  of  the  ocean  independently 
of  the  generating  mechanism  using  a  realistic  density  profile. 

The  model  used  in  this  study  was  derived  by  Rooth  (1970)  and  is 
presented  fully  in  Appendix  B.   It  makes  use  of  the  non-linear,  time 
dependent  equations  of  motion  in  deriving  the  vorticity  equation. 
Hydrostatic  equilibrium  is  assumed  for  the  formation  of  the  divergence 
equation.   It  is  assumed  that  the  data  were  collected  a  sufficient 
length  of  time  after  the  storm  passed  to  ensure  that  motion  was 
essentially  geostrophic.   Combining  the  divergence  equation  and  the 
vorticity  equation  results  in  a  governing  equation  expressing  the 
conservation  of  potential  vorticity  in  a  geostrophic  velocity  field, 

v2§  +    *L  <££  =  o.  n-i 

§  is  the  lifting  of  a  density  surface  from  its  undisturbed  state  (called 
the  matching  level),  f  is  the  Coriolis  parameter,  and  N^(Zq)  is  the 
square  of  the  Brunt-Vaisala  frequency  of  the  undisturbed  state.   Zq  is 
the  geometric  height  in  the  undisturbed  state.   The  observed  mass 
distribution  perturbed  from  the  undisturbed  state.   The  equation  is 
the  result  of  a  transformation  from  cartesian  coordinates  with  geometric 
height  as  the  independent  variable  in  the  vertical  to  a  "natural" 
coordinate  system  where  density  is  the  vertical  coordinate.   Figure  2 
illustrates  the  relationship  between  §  and  Zq.   Geometric  depth  is  the 
sum  of  Zq  and  \.      Formally,  Z  =  Zq(p)  +  §  (Zq,t)  ,   where  r  is  the  radius. 
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Figure    2.      The    relation  between   density    ( P)  ,    undisturbed   depth    (Zq)    and 
stretching    (£)  . 


Not  only  is  the  numerical  solution  considerably  simplified  in  this 
formulation,  but  solving  in  terms  of  the  stretching  is  intuitively 
appealing.   It  is  the  stretching  which  generates  the  relative  vorticity. 

The  horizontal  Laplacian  of  §  transforms  directly  from  cartesian 
to  cylindrical  coordinates.  Equation  II-l  is  therefore  the  governing 
equation  in  the  latter  system  also.   Let 

§(ZQ,r)  =  R(r)Z(Z0) 
in  order  to  separate  variables ,  and  II-l  becomes 

,  f2   ZMR  +  ZV2R  =  0.  II-2 

PTzo) 

2  2 

For  notational  convenience,  N  (Zn)  will  be  written  as  N~. 


Rearranging  II-2  yields 


f 2  •  VL   =  ~v2r 
nJ   Z     R 


II-3 


Since  a  function  of  Zq  can  be  equal  to  a  function  of  r  only  if  they  are 
each  constants,  then 


4  .  11     =   -V2R   =  k2  TI_4 

N|   Z      ^ 

Equation  II -4  implies 

lL  Z"   -   k2Z  =  0.  II-5 

Equation  II-5  must  be  solved  numerically  because  of  the  inclusion  of 
the  Brunt-Vaisala  frequency  as  a  function  of  Zq,   Note  that  the 
selection  of  sign  for  the  separation  constant  implies  that  solutions 
are  sought  which  are  exponential  in  form  in  the  vertical  and  harmonic 
in  the  horizontal.   Equation  II-4  implies 

V2R  +  k2R  =  o.  II -6 

Writing  the  Laplacian  in  full,  one  obtains 


R"  +  R'+  k2R  =  0,  II-7 

r 

assuming  axisymmetric  motion.   Equation  II-7  is  a  form  of  Bessel's 
equation  and  has  the  solution 

R  =  B  J0(kr)  +  C  NQ  (k2) .  II-8 

Nq  is  a  Neumann  function  or  Bessel  function  of  the  second  kind.   Since 
the  solution  must  be  finite  at  r  =0,  the  constant  C  equals  zero. 
With  no  loss  of  generality,  the  radial  solution  is  normalized  to  unit 
amplitude  at  r  =  0.   The  separation,  then,  is  the  following: 

00 

§  =n?l  50  (z0>°)  Jo  (knr>-  IX"9 

There  is  an  infinity  of  solutions  dependent  on  the  value  of  l^ ,  the 

arbitrary  constant  of  separation  called  the  wavenumber.   In  the  appli- 
cation of  equation  II-9  to  this  study,  the  solution  can  be  reduced  to 
the  response  to  a  single,  dominant  wavenumber.   The  technique  will  be 
discussed  following  a  summary  of  the  numerical  method  used  to  solve 
equation  II-l. 

Equation  II-l  was  solved  numerically  on  the  IBM  370/155  computer 
at  the  University  of  Miami.   A  pre-programmed  scientific  subroutine 
package  was  utilized  which  is  a  modification  by  Hamming  of  Milne's 
predictor-corrector  method.   Since  the  predictor-corrector  method  is 
not  self-starting,  a  fourth  order  Runge-Kutta  starting  procedure  is 
used  to  generate  sufficient  points  (Ralston,  1965).   The  numerical 
method  requires  specification  of  £  and  B^/SZq  as  initial  conditions. 
At  the  bottom  the  stretching  of  an  isopycnal  is  zero.   It  is  assumed 
that  a  dominant  wavenumber  for  the  disturbance  exists  in  the  interior. 
It  will  be  shown  later  that  this  assumption  is  realistic.   There  are, 
then,  two  arbitrary  parameters,  wavenumber  and  slope  at  the  bottom, 
which  have  to  be  determined.   A  solution  for  specific  values  of  slope 
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and  wavenumber  can  be  obtained  by  devising  a  scheme  which  matches  the 
numerical  solution  to  the  density  structure  in  two  dimensions.   The 
matching  uses  up  the  two  degrees  of  freedom  associated  with  the 
arbitrary  parameters.  A  numerical  scheme  could  be  devised  which  would 
optimize  the  matching  process  by  some  best  fit  method.   It  was  felt, 
however,  that  a  subjective  graphical  method  was  justified  considering 
the  quality  of  the  data.   The  approach  chosen  was  to  select  an  inter- 
mediate density  surface  and  vary  the  input  to  the  computer  program 
until  the  amplitude  of  the  response  and  the  decay  of  the  jQ(kr)  Bessel 
function  matched  the  observed  perturbation.   Essentially,  the  numerical 
solution,  generated  as  an  initial  value  problem,  was  forced  to  meet 
conditions  on  a  surface  on  which  density  is  conserved. 

The  governing  equation  (II-l)  requires  the  Brunt-Vaisala  frequency 
of  an  undisturbed  ocean  as  input.   The  input  was  taken  from  the 
temperature  and  salinity  tables  from  R/V  ATLANTIS  stations  5446  and 
5447  (Fuglister,  1960).   The  stations  were  selected  because  they  display 
temperature-salinity  characteristics  fairly  typical  of  the  Sargasso 
Sea  and  are  relatively  unperturbed  by  local  variations.   Figure  3  shows 
the  observed  response  of  the  ocean  to  the  storm. 
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Figure   3.      Observed  mass   structure   of   ocean  after   passage    of  Hurricane 
Kara.      Isopleths    are   values   of   sigma-t    (a    )  . 


III.   DISCUSSION  OF  RESULTS 


A.   DISCUSSION  OF  SOLUTION  AND 
CONSISTENCY  OF  RESULTS 


Figure  4  shows  a  comparison  between  the  computed  and  observed 
curves.   The  observed  mass  distribution  has  been  shifted  to  a  common 
vertical  axis.   This  step  is  consistent  with  the  assumption  of 
axisymmetric  motion  implied  in  the  transformation  of  the  Laplacian 
in  equation  II-6  to  cylindrical  coordinates.  As  seen  in  Figure  3,  the 
axis  of  the  observed  curves  slants  slightly  away  from  the  vertical. 
With  the  exception  of  the  near  surface  layer  and  a  few  isolated  spots 
at  the  sides,  where  the  assumption  of  axisymmetric  motion  is  in  doubt, 
the  difference  between  the  curves  is  no  more  than  10%  of  the  total 
perturbation  amplitude  (trough  to  crest)  of  any  isopycnal.   In  most 
cases  the  percentage  is  considerably  lower.   The  cross-hatched  area 
in  Figure  4  indicates  the  areas  where  the  deviation  of  the  two  curves 
is  significant.   Considerable  flattening  of  the  density  surfaces  occurs 
at  least  as  deep  as  230  meters ,  and  there  is  no  evidence  of  upwelling 
above  100  meters. 

The  lack  of  information  about  the  state  of  the  near  surface  layer 
prior  to  and  during  the  storm  precludes  a  complete  analysis  of  the 
causes  of  the  deviation  between  the  two  curves.  According  to  Fuglister 
(1947)  ,  average  surface  temperatures  for  October  in  the  area  of  the 
storm  range  from  25.8°C  to  26.1°C.   These  temperatures  correspond  very 
well  to  the  surface  temperatures  found  after  passage  of  Kara.   If  the 
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Figure  4.   Comparison  between  the  observed  mass  structure  and  the  computed 
response  after  Hurricane  Kara.   jQ(kr)  =  0  at  80  km. 
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near  surface  layer  did  not  cool  significantly,  then  one  must  ask  if 
upwelled  water  penetrated  the  thermocline  and  was  mixed  in  the  upper 
layer.   It  is  possible  that  the  apparent  lack  of  surface  cooling  in 
the  presence  of  obvious  upwelling  is  the  result  of  either  horizontal 
advection  in  the  surface  layer  moving  the  cold  water  away  from  the  place 
of  formation, or  advection  in  the  interior.   The  former  is  more  likely, 
but  the  latter  can  partially  explain  the  slanting  of  the  axis  of 
rotation  of  the  vortex  and  the  fact  that  the  dome  is  located  significantly 
to  the  south  of  the  published  hurricane  track. 

Upwelling  clearly  took  place  as  a  response  to  Kara.  One  would 
expect  to  find  evidence  of  downwelling  at  the  sides  to  satisfy  continuity 
requirements.   O'Brien  (1967)  finds  downwelling  in  a  numerical  model 
and  Leipper  (1967)  finds  evidence  of  it  from  the  changes  in  shape  and 
depth  of  the  thermocline.   Surface  data  taken  after  Kara  do  not  indicate 
downwelling.   The  mixed  layer  depth  is  remarkably  constant  across  the 
vortex  and  is  typical  of  the  ocean  in  this  region.  At  the  sides  of  the 
vortex,  in  the  interior,  downwelling  definitely  occurred  in  balance  of 
the  upwelling  in  the  center  region.   This  can  be  seen  by  the  dipping 
of  the  isopycnals  below  the  undisturbed  level.   Ths  disparity  between 
the  curves  at  the  top  of  the  pycnocline  suggests  significant  downward 
penetration  of  momentum  below  the  depth  of  the  mixed  layer,  leading  to 
an  effective  Ekman  depth  of  about  230  meters.   That  is,  if  the  Ekman 
layer  is  considered  that  region  where  friction  effects  dominate,  and 
where  the  conservation  equations  are  not  applicable,  then  the  depth 
of  significant  differences  between  the  observed  and  computed  curves 
may  be  identified  as  the  Ekman  depth. 


15 


It   is    possible   to   check  the  consistency  of   the  results  with 

existing  ideas    of  surface  wind  stress  by  working  backwards    from  the 

theoretical  displacement   of   the   isopycnals.      The   technique  begins  with 

the  Ekman  pumping  relation 

pfw  =  curl  t ,  III-l 

where  p  is  density,  w  is  vertical  velocity,  and  t  is  wind  stress. 

Stokes'  theorem  may  be  written  as 

Jt.  dS  ■  JJcurl  T-dA,  III-2 

S        A 

or,  by  substituting  equation  III-l  for  curl  T,  as 

JV-dS  =  pJJwf-dA.  Ill -3 

S        A 

Assuming  an  axisymmetric  wind  stress  distribution,  equation  III-3 

becomes 

pwfrrr2  =  2nrT  ,  III -4 

or 

t    =  wfr/2.  III-5 

If  w  =  6?/6t  where  6§  is  the  displacement  of  the  most  perturbed 

density  surface  and  6t  is  the  time  the  storm  was  stationary,  then  the 

average  stress  distribution  is 

t  =  6|-fr  .  III-6 

2-6t 

The  average  stress  may  be  found  by  taking  an  average  over  the  curve 

§  =  §0  Jo<kr) •  Then 

R  2tt 
t  =  I  •  I   f  f  g-r-d0-dr  m.y 

2n   R  0Q   2-6t 

.  f  .  |0  ?  f   JQ(kr)-r-d0-dr 
2tt   R  q  Jj       2-  6t 

To  simplify  the  integration,  it  was  decided  to  utilize  the  series  form 
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of  the  Bessel  function: 

Z2 


Noting  that 


and  us  ing 


J°(Z)=1"IW+  2W  +2*W+-"   •   IIr"9 


2tt 

U  d0  =  1, 
2tt  ft 


X  =  kr,  dr  =  dX/k 


to  change  variables  and  integration  limits ,  then  III-8  becomes 

k2-R  §   2'6t 
The  dominant  wavenumber  found  by  matching  the  theoretical  and  observed 
curves  is  3.0  x  10"->.   Therefore,  the  zero  crossing  occurs  at  r  =  80  km. 
Using  the  well  known  relation 

"  "  Pa%  =  Pwv*  =  CDU10  '  pa  "I"" 

enables  one  to  use  the  average  stress  calculated  in  equation  111-10  to 

find  wind  velocity,  or 

U?_  =   T  111-12 

10   CD-Pa 

where  u  is  the  friction  velocity  in  air,  v  is  the  friction  velocity 
in  water,  p   is  the  density  of  air  (1.3  x  10~J  gm  cm~3) ,  pw  is  the 
density  of  water  (1  gm  cm~3)  ,  CD  is  the  drag  coefficient,  and  U^q  is 
the  wind  velocity  at  a  specified  height.   The  wind  velocity  is  highly 
dependent  on  the  choice  of  drag  coefficient.   Unfortunately,  there  are 
numerous  estimates  of  the  dependence  of  drag  coefficients  on  wind  speed 
and  few  agree.   From  Roll  (1965)  there  seems  to  be  a  bunching  of 
estimates  for  near  hurricane  force  winds  at  C-q  =   2.5  x  10"^.   This  is  a 
figure  O'Brien  (1967)  also  uses  and  terms  "classical."  The  results  of 
this  calculation  for  various  values  of  the  drag  coefficient  are 
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tabulated  below.  C-q   =  1.5  x  10~3  is  a  widely  quoted  value  which  Malkus 
and  Riehl  (1960)  found  worked  well  in  their  model.   2.0  x  10-3  is  the 
mean  of  the  two  values. 

CD  U10(KTS) 


1.5  x   10"3 

57.6 

2.0  x   10"3 

50.2 

2.5  x   10"3 

45.2 

Air  Force  penetrations  of  the  storm  during  the  time  it  was  quasi- 
stationary  found  wind  speeds  between  45  and  55  knots.   According  to 
Colon  (1966),  the  wind  distribution  in  the  lower  few  thousand  feet  is 
uniform,  or  nearly  so.   Therefore  the  wind  speeds  should  be  repre- 
sentative of  those  near  the  surface.   The  tabulated  results  of  the 
calculations  are  inconclusive.   One  may  say  there  is  no  inconsistency 
between  the  theoretical  result  and  gathered  data.   Occasionally  the 
results  of  theoretical  models  are  used  to  fix  parameters  upon  which 
the  models  critically  depend.  A  definite  statement  cannot  be  made  in 
this  case.   However,  the  results  indicate  that  the  drag  coefficient 
probably  falls  within  the  above  range  for  winds  of  near  hurricane  force. 

As  will  be  shown  later  in  the  paper,  the  penetration  of  the 
response  of  high  wavenumber  motion  is  less  than  that  of  lower  wave- 
numbers  due  to  the  effect  of  the  wavenumber  on  the  exponential-like 
solution  in  the  vertical  (equation  II-5) .   The  ocean  acts  as  a  low- 
pass  filter  to  the  response.   This  means  that  the  difference  between 
the  two  curves  at  the  level  where  they  diverge  significantly  may  be 
due  partially  to  the  effect  of  the  response  at  higher  wavenumbers , 
which  will  change  the  height  of  a  density  surface  in  some  unknown  fashion 
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It  is  also  possible  that  the  presumed  center  of  the  vortex  was 
not  the  exact  center.   If  not,  then  missing  the  center  can  account  for 
some  of  the  difference  between  the  curves. 


B.   THE  MIXING  PROCESS 
AND  MOMENTUM  TRANSFER 


Several  calculations  can  be  made  which  place  an  upper  limit  on 
the  effects  of  the  mixing  process.   The  calculations  assume  a  balance 
in  the  vertical  between  the  downward  diffusion  of  heat  due  to  mechanical 
mixing  and  the  upwelling  of  cold  water  from  below  the  local  thermo- 
cline.   The  governing  equation  for  this  process  is 

w-St/SZ  =<-o2T/az2,  111-13 

where  w  is  vertical  velocity,  T  is  temperature  and  <  is  the  eddy 
diffusivity  of  heat.   The  solution  of  111-13  is 

T  =  T0  exp(w/<)Z.  111-14 

</w  is  the  e-folding  distance  for  the  temperature  decay  and  is  therefore 
the  scale  height  of  the  local  thermocline  in  this  problem.   H^,  the 
thermocline  scale  height,  may  be  found  by  taking  the  logarithm  of  111-14 
obtaining  the  following: 

&ff  =  2jiTq   +Z/Ht.  111-15 

Plotting  8/riZ   vs  Z  and  finding  the  slope  determines  the  scale  height. 
Equation  111-15  is  the  equation  of  a  straight  line  on  semi-log  paper, 
and  1/H-j  is  the  slope  of  the  line.   Figure  5  is  such  a  plot.   Since 
pure  mixing  is  assumed,  the  vertical  velocity  may  be  approximated  by 
6h/6t,  the  difference  in  lifting  between  the  observed  and  calculated 
curves  divided  by  the  relevant  time  scale  of  one  day.  <   determined  in 
this  manner  is  1400  cm  /sec.   This  value  of  <  is  several  orders  of 
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Figure   5.      Log-linear   plot   of   temperature    (°C)    vs    depth    (m) 
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magnitude  too  large  for  abyssal  processes  but  is  probably  the  correct 
order  for  surface  processes  where  eddy  viscosity  and  eddy  diffusivity 
may  be  nearly  the  same.   If  the  point  of  divergence  of  the  curves  is 
at  230  meters,  the  eddy  viscosity  found  from  the  Ekman  depth  formula, 
do  =  Tf\j2\)/f    ,  is  2000  cm  /sec.   1400  cm2/sec  is  almost  certainly  an 
overestimate  of  the  eddy  diffusivity  since  Ekman  divergence  effects 
are  neglected.   If  the  depth  of  the  mixed  layer  is  neglected  in  the 
calculation  of  the  eddy  viscosity,  dg  =  180  meters  and  v  =  1200  cm2/sec. 
This  calculation  assumes  the  existence  of  two  separate  regimes:   the 
mixed  layer  where  eddy  viscosity  coefficients  are  very  high;  and  the 
layer  of  Ekman  divergence  below  the  frictional  layer,  indicated  by  the 
disparity  between  observed  and  computed  results.   The  similarity  between 
the  values  of  eddy  viscosity  and  eddy  diffusivity  is  presented  for 
completeness  and  to  suggest  topics  for  further  investigation. 

Under  conditions  of  negligible  advection,  it  is  possible  to 
extend  an  energy  calculation  by  Turner  (1969).  Assuming  pure  mixing, 
the  change  in  potential  energy  between  the  observed  and  theoretical 
curves  can  be  used  to  work  backwards  to  find  wind  velocity.   Selecting 
an  isopycnal  where  the  two  curves  do  not  diverge  and  calling  the  density 
Pq  a  reference  density,  one  can  plot  Pz~Po  vs  Z2  (depth  squared). 
Taking  one  half  the  area  under  the  resulting  curve  solves  the  energy 
equation, 

EP0T  =  gJ(Pz-Po)Z  d  Z'  111-16 

graphically.   Figure  6  is  a  plot  of  Ap  vs  Z^  for  the  center  station, 
where  maximum  upwelling  occurred.   The  change  in  potential  energy  found 
planimetrically  is  10.14  x  107  ergs  cm"2.   By  dimensional  arguments, 
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Figure  6.   Density  difference  (gm/cm3)  plotted  vs  square  of  depth  (cm2) 
for  potential  energy  calculation.   Station  #24  is  compared  to  the  center 
of  the  computed  curves . 
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Turner  asserts  that  it  is  reasonable  to  write  the  time  rate  of  change 
of  energy  as 

dEp0T  -  C^Uxq  =  C2PU103.  111-17 

dt 

The  tabulation  below  gives  values  of  wind  speed  for  two  values  of  the 

constant  Co  found  by  Turner. 

C2 Uiq(KTS) 

2.6  x  10"8        71.4 

1.3  x  10"8        90 


Both  values  of  wind  speed  are  too  high  if  pure  mixing  is  responsible 
for  the  potential  energy  change.   If  one  accepts  50  knots  as  a 
reasonable  value  of  wind  speed  during  the  storm,  then  surface  mixing 
accounts  for  about  half  the  energy  change.   That  is,  at  least  half  of 
the  density  field  discrepancy  from  the  theoretical  curves  must  derive 
from  the  deep  Ekman  divergence.  As  the  drag  coefficient,  Turner's 
constants  are  not  well  defined.   The  dependence  of  Turner's  constants 
and  the  drag  coefficient  on  wind  speed  is  unknown.   Turner's  method 
requires  a  one -dimensional  process  and  a  regime  in  which  surface  cooling 
and  resultant  convective  overturning  are  unimportant.   Both  processes 
are  usually  significant  under  hurricane  conditions,  although  their 
effect  may  be  minimal  in  this  case  since  surface  cooling  apparently  did 
not  occur.   Selecting  the  center  station  minimizes  the  advective  effects, 
since  §  is  a  maximum  there  and  d§/5r  approaches  zero.  As  Turner 
correctly  points  out,  his  method  is  illustrative  but  hardly  conclusive. 

C.   THE  EFFECT  OF  LENGTH  SCALE  ON  RESPONSE 

Some  interesting  properties  of  the  solution  as  a  function  of  scale 
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of  forcing  and  depth  are  illustrated  in  Figure  7.   The  governing  equation 
(II-l)  was  solved  using  the  same  initial  conditions  and  density  profile 
but  different  values  of  wavenumber.   The  Length  scale  was  determined 
by  selecting  the  value  of  radius  which  corresponded  to  the  first  zero 
of  the  Jq  Bessel  function,  found  from  the  formula  (1 /wavenumber)  x  2.4, 
The  values  of  §  were  normalized  to  100  meters  of  Ekman  pumping  at 
Zq  =  25  meters  and  isopleths  of  constant  §  were  plotted.   The  response 
may  be  split  into  a  barotropic  and  baroclinic  regime.   In  the  barotropic 
regime,  response  is  characterized  by  the  stretching  being  a  linear 
function  of  depth,  since  stretching  is  proportional  to  vertical  velocity. 

The  transition  from  a  baroclinic  to  a  barotropic  response  is  rather 
broad  in  scale,  occurring  between  300  to  500  km.   For  scales  greater 
than  500  km,  the  response  may  be  continued  to  be  barotropic.   The 
response  for  high  wave  numbers  is  confined  to  a  shallow  surface  layer. 
The  latter  provides  the  rationale  for  basing  the  solution  of  the 
vorticity  equation  on  the  existence  of  a  dominant  wavenumber  in  the 
interior.   The  dependence  of  the  vertical  structure  of  the  response  on 
the  scale  of  the  motion  suggests  a  definite  problem  in  the  representation 
of  a  continuously  stratified  ocean  by  a  two-layer  model.   For  the  higher 
wavenumbers ,  the  two  layer  model  lacks  sufficient  resolution  to  describe 
the  response  adequately. 

The  ordinate  in  Figure  7  is  undisturbed  depth  (Zq) •   Since  depth 
is  equal  to  the  sum  of  the  undisturbed  depth  and  the  stretching,  zero 
depth  in  the  figure  should  not  be  confused  with  the  free  surface.   The 
free  surface  lies  considerably  below  Zq  =  0. 

It  is  anticipated  that  similar  diagrams  can  be  generated  using 
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Figure  7.   Isopleths  of  stretching,  £ (m)  ,  as  a  function  of  depth  and 
length  scale.   Values  of  E,   are  normalized  to  100  meters  pumping  at 
25  meters  matching  depth. 
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density  distributions  typical  of  other  areas  of  the  ocean.   The  diagrams 
should  prove  useful  in  predicting  the  response  in  these  areas .  A 
density  profile  with  higher  resolution  than  that  obtained  from 
Fuglister's  Atlantic  Ocean  Atlas  is  preferable. 

D.   GEOS TROPHIC  VELOCITIES 

It  is  possible  to  calculate  geostrophic  velocities  using  the 
solution  of  equation  II-l.   Noting 

Z(p,r)  =  Z0(p)  +?(ZQ,r),  111-20 

and    from  Appendix  A, 

9p\        =-3z\     dp_ 
3rjZo   "   drjp    SZ0      ' 
one  may  use   the    thermal  wind   equation, 


111-21 


M   =  z&.  m\       ,  in-22 

3z0     pf  aryz0 

to  find  velocities  throughout  the  column  and  for  any  radius. 

Substituting  111-21  into  111-22,  one  obtains 

|U   =  _g.  |z\   Jg.  .  111-23 

dZ0   Pf  °*JP  ^0 

Equation  111-23  may  be  integrated  in  the  vertical: 

j^du-i;0*!*  -(%\  dz0. 

uz^  f  zy  p  ^o     V3rJp      ° 

Noting  Uzn    (bottom)    =  0  and  Nj*  -  Z&  M  ,   equation  III -24  becomes 
0  u         p    oZq 

Introduce   111-20  for  Z    in  111-25,   obtaining 


111-24 


f  £"     u     \Sr/p 
^0 

Geostrophic  velocities  were   calculated  numerically  using  equation  111-26 
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and  the  result  established  previously,  §  =  §q  jQ(kr) .   The  derivative 
of  §  is 

d^/Br  =  k  §0  J1(kr) .  111-27 

Equation  111-26  is  a  valid  expression  for  the  velocities  induced  by 
Kara,  except  in  the  Ekman  layer,  where  the  predicted  density  field 
departs  from  the  observed.  Veronis  (1956)  showed  that  motion  which 
has  forcing  on  scales  larger  than  the  baroclinic  radius  of  deformation 
is  essentially  geostrophic.   In  this  case,  the  scale  of  forcing,  80  km, 
is  more  than  twice  the  baroclinic  radius  of  deformation,  36  km. 
Figure  8  is  a  comparison  of  the  geostrophic  velocities  based  on 
equation  111-26  and  velocities  calculated  using  the  dynamic  method 
(called  the  observed  velocities) .   The  reference  level  for  calculating 
the  observed  velocities  was  found  using  Defant's  method.   The  level 
was  conspicuously  present  in  all  the  profiles.  Although  the  percentage 
of  error  is  high  at  the  lower  depths,  the  absolute  magnitude  difference 
is  low.   In  general,  the  magnitude  of  error  remains  the  same  as  depth 
decreases.   The  two  velocity  profiles  agree  rather  well  in  magnitude 
and  general  shape.   Aside  from  the  reasons  for  the  difference  between 
the  curves  in  the  Ekman  layer,  which  have  been  discussed,  the  calculation 
of  geostrophic  velocities  poses  special  problems.   In  this  case,  the 
distances  between  stations  used  to  calculate  the  velocities  by  the 
dynamic  method  were  extreme.   It  was  found  that  velocities  at  a  given 
radius  and  depth  varied  as  much  as  5  cm/sec  depending  on  the  particular 
combination  of  stations  used.   The  curves  represent  the  smoothest 
combination.   The  observed  velocities  were  computed  from  the  non-centered 
observed  data,  which  certainly  accounts  for  some  of  the  difference 
between  the  curves.   The  effect  of  centripetal  acceleration  has  been 
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Figure  8.   Comparison  of  numerically  calculated  geostrophic  velocities 
(cm/sec)  with  velocities  found  using  the  dynamic  method  (observed) . 
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neglected  in  the  calculation  of  velocities,  and  in  the  formulation  of 

the  vorticity  conservation  model.   The  inclusion  of  centripetal 

accelerations  makes  the  model  non-linear,  introducing  an  undesired 

degree  of  difficulty  in  solving  the  vorticity  equation.  A  measure  of 

the  importance  of  the  centrifugal  force  term  is  the  Rossby  number,  v 

f -r 

(Neumann  and  Pierson,  1966)  .   Choosing  r  =  40  km,  v  =  40  cm/sec,  and 
f  =  7.5  x  10~5  sec"!,  the  Rossby  number  is  0.13.   It  can  be  shown  that 
geos trophic  velocity  is  an  overestimate  of  the  gradient  velocity  if 
curvature  is  important  (Haltiner  and  Martin,  1957). 

The  success  in  obtaining  substantial  agreement  between  the 
theoretical  mass  field  and  velocities  and  the  observations  supports  the 
validity  of  the  initial  model  formulation. 


IV.   SUMMARY  AND  CONCLUSIONS 

A  model  based  on  the  conservation  of  potential  vorticity  was 
solved  numerically.   The  solution  was  utilized  to  investigate  some 
aspects  of  the  oceanic  response  to  a  tropical  storm  as  a  function  of 
horizontal  scale.   In  particular,  the  interior  response  of  a  moderately 
intense,  quasi-stationary  cyclone  was  duplicated.   The  matching  was 
successful  below  a  depth  of  230  meters  which  was  associated  with  the 
Ekman  depth.   The  lack  of  information  concerning  the  state  of  the 
ocean  in  the  near-surface  layer  prior  to  the  passage  of  Hurricane  Kara 
compounded  the  problem  of  discovering  what  caused  the  divergence  of 
the  observed  and  computed  curves.   The  divergence  was  attributed  to  the 
neglect  of  high  wavenumber  modes  when  solving  for  the  interior  response  , 
and  to  some  ill-defined  combination  of  vertical  mixing  and  momentum 
transfer  by  friction.   Gross  estimates  indicate  that  the  two  processes 
are  of  equal  importance  in  this  particular  case.  An  upper  limit  for 
the  eddy  diffusivity  is  1400  cm  /sec.   Vigorous  mixing,  if  it  occurred, 
may  not  have  resulted  in  significant  cooling  of  the  mixed  layer. 

Geostrophic  velocities  were  calculated  from  the  raw  data  using 
Def ant's  method  to  find  a  reference  level  of  no  motion.   These 
velocities  were  compared  to  those  calculated  from  the  numerical  solution. 
The  velocities  agree  rather  well  below  the  near -surface  layer.   Since 
the  velocity  computation  integrates  the  stretching,  the  comparison  of 
velocities  should  not  be  expected  to  be  as  successful  as  the  comparison 
of  §  itself.   The  integration  of  §  means  that  errors  are  accumulated 
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over  the  entire  column.  Geostrophic  vertical  shear  matches  better 
than  geostrophic  velocities.   It  is  plausible  that  the  nearly  constant 
error  in  magnitude  between  the  velocity  fields  is  a  result  of 
inaccuracies  in  the  dynamic  computations;  specifically,  limitations  in 
the  Defant  method  of  obtaining  the  level  of  no  motion  may  be  the  cause. 

It  was  found  that  the  scale  of  forcing  can  have  dramatic  effects 
on  the  penetration  of  the  response.   For  equal  Ekman  pumping,  higher 
wavenumber  (short  scale)  response  does  not  penetrate  very  deep  into 
the  ocean  compared  to  low  wavenumber  response.   Thus  the  scale  of  the 
disturbance  must  be  considered  when  choosing  the  depth  and  number  of 
layers  for  modelling  a  continuously  stratified  fluid.   A  two  layer 
model  may  lack  the  resolution  necessary  to  accurately  describe  the 
response  in  the  upper  layer.  Adequate  modelling  of  the  near-surface 
layer  by  a  model  such  as  the  one  used  in  this  study  requires  the 
inclusion  of  high  wavenumber  motion. 
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APPENDIX  A 


THE  NATURAL  COORDINATE  SYSTEM 
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Although  the  standard  cartesian  coordinate  system  is  the  one 
normally  used  in  hydrodynamics ,  there  is  no  reason  why  the  vertical 
coordinate  must  be  geometric  height,  a  fact  recognized  long  ago  by 
meteorologists  and  some  oceanographers .   The  use  of  a  "natural"  vertical 
coordinate,  such  as  pressure,  density  or  potential  density,  can  simplify 
an  analysis  but  render  boundary  conditions  formidable. 

In  the  potential  vorticity  conservation  model,  density  will  be 
utilized  as  the  "natural"  vertical  coordinate.   This  is  advantageous 
because  it  enables  one  to  solve  directly  for  the  displacement  of  an 
isopycnal  from  equilibrium.   The  method  of  coordinate  transformation  is 
discussed  in  Haltiner  and  Martin  (1957) .   It  is  summarized  below  since 
the  method  is  implied  in  the  model  and  is  used  directly  in  deriving  the 
geostrophic  velocities  (equation  111-26) . 

Consider  0(x,y,z,t)  =  constant.   The  change  of  Tl(x,y,z,t)  along 

0  =  constant  is  equal  to  the  change  along  the  horizontal  plus  the 

change  along  the  vertical  due  to  the  slope  of  0,  i.e.   From  Figure  9, 

6lA    =6^    +  |Tl\   .  A-l 

/a-c     Ja-b  p-c 


Now, 


/y>0>to 


=  |lh  6x    +  JJfl  6z  A-2 

dx/y,z,t0   oz/x,y,t0 


by  the  chain  rule  at  a  constant  time,  tQ.     A-2  may  be  written  as 

M\         -  jnj\         +   v$\  2sN  ,         A-3 

dx/y,0,to     oX//y,z,t0       dzyx,y,t0  dx/y,0,to 

where  ^2.  is    the   slope  of  the  0  surface   in   the  x  direction.      Let  T|  =  p 
ox 

and  A-3  becomes 

*b\  =  2eN   .  *£.  .  h.  a-4 

3xyz   9x/0   3z   3x 
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Figure  9.   0  =  constant  surface  separated  into  horizontal  and  vertical 
components  to  illustrate  coordinate  transformation  technique. 
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APPENDIX  B 


THE  POTENTIAL  VORTICITY  CONSERVATION  MODEL 
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The  governing  equation  for  the  interior  motion  of  a  hurricane,  it 
is  hypothesized,  expresses  the  conservation  of  potential  vorticity. 
The  equations  of  motion  used  are  non-linear,  time  dependent,  hydrostatic 
equations  in  a  rotating  system  characterized  by  the  Coriolis  parameter, 
f,  which  is  a  function  of  the  horizontal  coordinates  alone: 

B-l 


u,.  +  u  •  Vu  -  fv  +  TT   =  0 


V.  +  U  *  Vv  +  fu  +  TT   =  0, 


B-2 


where  tt  =  *-  +  gZ  .   The  vorticity  equation  may  be  formed  by  cross 
differentiating  equations  B-l  and  B-2,  resulting  in  B-3: 


a         a    j.     a 

_  +  U  4j—  +  v  ^~ 

dt      "X      Oy 
B-3  may  be  written  as 

3__  +   3  +  v  a_ 

3t     ox     Oy 


(vx-uy)+(f-»vx-uy)Cux-tvy)+  ufx  +  vfy  =  0  B-3 


(f\-0+(^v-"J(Uv\)   ■  0» 


x  y 


since  f  is  not  a  function  of  time.   That  is, 


d 
dt 


(f-tvx-uy)(ux+vy)  =  0, 


B-4 


B-5 


where 


dt  —  dt     Ox     Oy 
The  divergence  equation  may  be  formed  by  taking  r—  (equation  B-4) 


Bx 


and  adding  to  _  (equation  B-2) . 

3y 


.2.2. 


§j    (ux4vy)-f(vx-uy)    +  V^TT  +  fyu-fxvhax^4vxuy^yvx    =   0        B-6 
Assuming  geostrophic  motion  reduces    the   balance  equation  to 


f(v   -u  )    =  V2tt. 

x     X        y' 

The  continuity  equation  is 


u„+v  +w„  =   0. 

X   y   Z 


B-7 


B-8 


In  the  vertical,  w  depends  only  on  z(p),  therefore 
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(u+v   )  +  Hw  .  *£  =  0. 

x    y        op   oz 

Noting  w  =  dz/dt,  equation  B-8  becomes 

op   x  y     dt\oP; 

Substitute  equation  B-9  into  the  vorticity  equation  (B-4)  to  obtain 
d  J-  (°z/5p) 

a-(f  +  v    -  u  )  -  (f  +  vx  -  uv>  •  H 1  =  o.  B-io 

dt       x     y  x     y       oz/op 

Rearranging   terms,   equation  B-10  becomes 

ir<f  +  v*  -  V       _  ^(°2/aP)  =  n  b_u 

f  +  VX  -  ^  ^Z/^P 

or,  by  definition  of  the  logarithm, 

X— * — >=  0.  B-12 


dt  \  dz/dp 

Equation  B-12  implies  that  for  each  fluid  element 

"^""y"1^  =  constant.  B-13 

oz/dp 

The  existence  of  a  uniform  reference  state  implies  that  the  constant  is 

a  function  of  density  alone.   If  the  reference  state  is  rest,  relative 

vorticity  is  zero,  and  the  constant,  0(p)  ,  equals  f/(3z/dp)Q.   The 

vorticity  equation,  then,  becomes 

^-V*  =  _f .  B-14 

3z/op   (oz/op)0 


Substituting  equation  B-7  into  B-12  and  rearranging,  one  obtains 

B-15 


l+V2nVdz\       -  dz 


i  /Wo 


f2  Aapyo     sp 

Next,  apply  the  hydrostatic  equation, 

TT  =  £■  +  gz.  B-16 

Rearranging  and  taking  the  partial  derivative  with  respect  to  z  results 


in 


|^(rr-gz)p  =PZ.  B-17 
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Assuming  hydrostatic  equilibrium,  B-17  becomes 

— (n-gz)p   =  -gp.  B_18 

oz 

Integrating  B-18  results  in  the  following: 

Z 
pTT-pgz  =  pQ  TT0  -  J  pgdz.  B-19 

0 

The  last  term  in  equation  B-19  can  be  integrated  by  parts: 
Z  p 

-Jpgdz  =  -pgz  +  Jgzdp  .  B-20 

0  PO 

Substituting  B-20  into  B-19,  the  desired  forms  of  the  hydrostatic 

equation  is  achieved: 

tt  =  tt0  +  I  Jgzdp  .  B-21 

P  PO 

Take  the  Laplacian  of  B-21: 

p 
V2tt  =  V2TTQ  +  I  JgV2zdp    .  B-22 

P  Po 

Substitute  equation  B-22  into  equation  B-15: 

aP   WJo     V      o     PJ«      w^Q 

Therefore,  if  Pq  is  chosen  as  the  surface  density,  where  h  is  the  free 
surface  height, 

32.  YisN     =  %(y2h  +  I  fv2zdp)(|5.V    .  B-23 

op  \OP;o    f2        p  JPo      W° 

The  depth,  z,   may  be  split  into  two  parts, 

Z    =  Ztfp)    +  §(x,y,z0). 

Equation  B-23  becomes,   then,   after  rearranging, 

V2h  +ljv2?dp    =£|i/k\      .  B-24 

P    p0  §      9P    W° 


Next,   take   the  derivative   of  B-24  with  respect   to  p: 

p  gpo  WV  z/°       gpo  Bp  \av 


B-25 
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Recall 

i-  =  (*£\    i_ 

dp  ""  V^P/O  3z0 

so  equation  B-25  becomes 

I  V2?  =  -f2„,    £l 
PO       PON^(Zo)  dZ02  ' 

or 

V2§+-^ l^r-  "  0.  B-26 

NZ(Z0)  &z02 

Equation  B-26  is  the  governing  equation  used  in  the  analysis  described 
in  the  text.  The  basic  derivation  is  attributed  to  Rooth  (1970). 
Several  steps  were  expanded  by  this  author  for  clarity  and  completeness. 
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